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ABSTRACT 

The nonlinear evolution of the Kelvin-Helmholtz instability is a popular test for code verifica¬ 
tion. To date, most Kelvin-Helmholtz problems discussed in the literature are ill-posed: they 
do not converge to any single solution with increasing resolution. This precludes comparisons 
among different codes and severely limits the utility of the Kelvin-Helmholtz instability as a 
test problem. The lack of a reference solution has led various authors to assert the accuracy 
of their simulations based on ad-hoc proxies, e.g., the existence of small-scale structures. 
This paper proposes well-posed Kelvin-Helmholtz problems with smooth initial conditions 
and explicit diffusion. We show that in many cases numerical errors/noise can seed spurious 
small-scale structure in Kelvin-Helmholtz problems. We demonstrate convergence to a ref¬ 
erence solution using both Athena, a Godunov code, and Dedalus, a pseudo-spectral code. 
Problems with constant initial density throughout the domain are relatively straightforward 
for both codes. However, problems with an initial density jump (which are the norm in 
astrophysical systems) exhibit rich behavior and are more computationally challenging. In 
the latter case, Athena simulations are prone to an instability of the inner rolled-up vortex; 
this instability is seeded by grid-scale errors introduced by the algorithm, and disappears as 
resolution increases. Both Athena and Dedalus exhibit late-time chaos. Inviscid simulations 
are riddled with extremely vigorous secondary instabilities which induce more mixing than 
simulations with explicit diffusion. Our results highlight the importance of running well-posed 
test problems with demonstrated convergence to a reference solution. To facilitate future 
comparisons, we include the resolved, converged solutions to the Kelvin-Helmholtz problems 
in this paper in machine-readable form. 
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1 INTRODUCTION 

The Kelvin-Helmholtz (KH) instability results from a wide array 
of velocity-shear profiles in a continuous fluid, or across the in¬ 
terface between two distinct fluids. The instability is ubiquitous 
in nature, playing important roles in meteorology, oceanography, 
and engineering. The KH instability plays a particularly prominent 
role in astrophysical systems ranging in scale from stellar interiors 
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(e. g. Bruggen & Hillebrandt 2001) and protoplanetary disks (e. g. 
Johansen et al. 2006) to the evolution of the intergalactic medium 
(e. g. Nulsen 1982, 1986). Physically, the KH instability wraps up 
coherent sheets of vorticity into smaller, less organized structures. 
The small scale motion then stretches and cascades to yet smaller 
scales. The instability therefore plays fundamental roles in fluid 
mixing and in the transition to turbulence. 

Because of its prevalence in nature and its physical signifi¬ 
cance, KH test problems are commonly used to evaluate the accu¬ 
racy of different astrophysical hydrodynamics codes (e. g. Springel 
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2010; Hopkins 2015; Schaal et al. 2015): if a code can properly 
simulate the KH instability, it is presumed to capture mixing and 
turbulence in astrophysical simulations. Ideally, such an important 
test problem should stand against an analytic solution to ensure the 
veracity (not just reproducibility) of simulation results. Some ana¬ 
lytic work addresses the KH instability with a sheet vortex model 
(Moore 1979), but only for incompressible fluid equations. For the 
compressible Navier-Stokes equations relevant to astrophysics, no 
analytic description of the nonlinear KH instability currently exists. 


Absent a nonlinear analytic prediction, a resolved reference 
simulation provides the only reasonable approximation of the true 
solution. Comparing to a well-controlled and high-resolution bench¬ 
mark gives a proxy for the true error of a given test. Robertson et al. 
(2010) and McNally et al. (2012) present careful studies of the 
early evolution of the KH instability. These authors also point out 
the numeric ill-posededness of contact-discontinuity simulations, 
in spite of existing analytical solutions in the linear and/or incom¬ 
pressible regimes. These works emphasize that converged nonlinear 
simulations require well-resolved initial conditions. One limitation 
of these studies, however, is that Robertson et al. (2010) and Mc¬ 
Nally et al. (2012) only provide converged reference simulations 
for the linear (and possibly weakly nonlinear) phase. In addition, 
converged nonlinear solutions require solving dissipative equations. 
Many available astrophysical codes do not implement this essential 
feature. As a result, these works could only follow the instability 
for a few e-folding timescales. 


Not all works take the benchmark approach, however. In place 
of a nonlinear reference solution, some authors use apparent small- 
scale structure as a proxy for the accuracy of their simulations 
(e.g., Springel 2010; Hopkins 2015). Presumedly, more small-scale 
structure implies less numerical dissipation, and therefore greater 
accuracy. We find in the current paper that this intuition can, in 
some cases, lead to false conclusions. Some tests also abandon the 
smooth initial conditions of Robertson et al. (2010) and McNally 
et al. (2012), even though this choice precludes convergence of even 
the linear phase of the instability because the linear growth rates 
increase with wavenumber for an initially discontinuous velocity 
profile. 


In this paper, we extend the work of McNally et al. (2012) 
by providing reference solutions for the strongly nonlinear evolu¬ 
tion of the KH instability. We use a smooth initial condition and 
explicit diffusion. We conduct simulations using both Athena (a 
Godunov code), and Dedalus (a pseudo-spectral code that can solve 
the Navier-Stokes equations of compressible hydrodynamics) and 
find that both converge to the same reference solutions. We see 
agreement among different codes and different resolutions, with 
the validity of the reference solution limited only by (unavoidable) 
chaotic evolution at late times. We propose that future code tests 
include this KH instability problem and compare to our validated, 
converged, reference solutions. 


We organize the remainder of the paper as follows. Section 2 
describes the equations, initial conditions, and codes used for our 
simulations. The results comprise two sections. In section 3.1 we 
discuss the simpler simulations with constant initial density. Sec¬ 
tion 3.2 discusses the more complicated simulations with an initial 
density jump. Section 4, summarizes our results. 


2 METHODS 

2.1 Equations and Initial Conditions 

We solve the hydrodynamic equations, including explicit terms for 
the diffusion of momentum and temperature: 

^+V-(p«) = 0, (la) 

ot 

^{pu)+ V • {P\ + pu ®u) --V'U, (lb) 

ot 

^ + V .[(£• + />)«]= V . (xpVT) - V . (« . □), (Ic) 

Ot 

along with the nondimensionalized ideal gas equation of state, P - 
pT, with constant ratio of specific heats y = 5/3. I is the identity 
tensor, is the thermal diffusivity (with units cm^/s; K - nk\^x is 
the thermal conductivity), and 

n = -vp -h {VuY - ^IV • wj (2) 


is the viscous stress tensor with viscosity v (with units cm^/s). We 
assume both y and;^ are constant. 

We add a passive scalar to our simulations which we refer to 
as “dye.” The local fraction of dye particles c expresses dye con¬ 
centration, and initially ranges from 0 to 1. The local conservation 
of dye is then 


d , ^ , X dc „ ^ 

— (pc) -h V • ipcu) = p— = -V • (g, 


dye’ 


Sdye ~ 


( 3 ) 

( 4 ) 


where djdt represents the Lagrangian derivative, and Vdye represents 
a diffusion coefficient for dye molecules (with units cm^/s). These 
equations conserve the total dye mass f pcdV. 

We define a dye entropy per unit mass s = - c Inc, along with 
its volume integral 


5 = 



sdV. 


These evolve such that: 



v-[(i+inc)e,y; 


dt 


IVcp 

- pydye - 

C 

r ivcP 

= J PVdye- dV > 0. 


( 5 ) 

( 6 ) 
( 7 ) 


The second term on the left-hand side of equation 6 represents 
the entropy flux due to reversible diffusion of the dye. The right- 
hand side represents entropy generation due to non-reversible dis¬ 
sipation.^ The volume-integrated entropy S satisfies the following 
important properties: 


(i) A fully unmixed fluid with c = 0 or c = 1 everywhere has 
zero entropy (S =0). 

(ii) A fully mixed fluid with c* = J pc dV/ f p dV maximizes 
the entropy, 5 max = -c* In c* Jp dV. 

(iii) S increases monotonically with time if Vdye > 0, and stays 
constant otherwise. 


We restrict our attention to periodic simulations. This avoids 


^ Equation 6 can be made to look like the analogous equation for heat 
conduction with the definition of a new “temperature” r^ye = - 
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potential difficulties with imposing Dirichlet and/or Neumann 
boundary conditions. Our initial conditions are: 


1 Ap 1 
p - \ + — X - I tanh 
Po 2 


= Wflow X |tanh ^j - tanh^-—- l| 


Uz- A sin(2;rx) X 


exp 




- exp 


{z - Zif 


P = Po 


c = - I tanh 




(8a) 

(8b) 

(8c) 

(8d) 

(8e) 


where a - 0.05 and cr - 0.2 are chosen so that the initial condition 
is resolved in all of our simulations. We take waow = 1 and Pq = 10 
so that the flow is subsonic with a Mach number M ~ 0.25. The 
size of the initial vertical velocity perturbation is A = 0.01. The 
Athena simulations are initialized with these functions evaluated 
at cell-centers even though Athena data represents cell-averaged 
quantities (see Appendix A for more discussion of this effect). 

We adopt a rectangular domain with x in [0, L^), and z in 
[0,2L^), with Ljc = I and 2L^ = 2, and zi = 0.5, Z 2 = 1-5, with 
periodic boundary conditions in both directions. The simulations 
have a horizontal resolution of N grid points (in Athena) or modes 
(in Dedalus) in the v direction, and 2N grid points/modes in the z 
direction. Our initial condition has a reflect-and-shift symmetry: 
taking z ^ 2 - z and v ^ x + 1/2 changes the sign of but 
leaves the other quantities invariant. Thus, the simulations solve 
for the same flow twice. This is a requirement when using periodic 
boundary conditions, but also provides a test of whether or not 
the numerical simulations can preserve the symmetry. Almost all 
simulations presented here maintain the symmetry. We therefore 
only show the lower half of the domain. We calculate volume- 
averaged quantities like the dye entropy or the L 2 norm with respect 
to the entire domain. 

In equation 8a, the free parameter Ap/po represents the density 
jump across the interface. We study simulations with Ap/po = 0 
in section 3.1 and with Ap/po = 1 in section 3.2. We refer to this 
change in density as a “jump” throughout, although the transition 
is smooth, set by the tanh in equation 8a. The Reynolds number Re 
quantifies diffusion. 


y=X = '^dye 


LAu 


(9) 


where Au = 2wflow is the change in velocity. Note that we set the 
thermal diffusivity x constant; consequently, the thermal conduc¬ 
tivity K oc p. Throughout the paper we measure time in units of 
T/Wflow. so t = 1 corresponds to approximately one turnover time. 
Equations 1-9 specify our system, with the free parameters Ap/po 
and Re. In the following section we detail our methods for solving 
this system of equations. 


2.2 Numerical Methods 

We study the KH instability using two open-source codes employ¬ 
ing very different numerical methods: Athena & Dedalus. 

Athena^ is a finite-volume Godunov code (Gardiner & Stone 
2008; Stone et al. 2008). The scheme represents all field quantities 
with volume averaged values in each grid element. A Riemann 
problem solves for fluxes between elements. We use third-order 

^ Athena is available at https: //trac. princeton. edu/Athena/. 


reconstruction with limiting in the characteristic variables to approx¬ 
imate field values at the element walls, the HLLC Riemann solver, 
and the CTU integrator. We used the “-03” compiler flag using In¬ 
tel 14.0.1.106 and Mvapich2 2.0b on the Stampede supercomputer. 
We repeated some runs using second-order reconstruction and/or 
the Roe Riemann solver and/or stricter compiler flags (e.g., “-02 
-fp-model strict”) — these choices did not qualitatively affect the 
solutions. We use a static, uniform mesh, and a CEL safety factor 
of 0.8. 

Athena is second-order accurate in both space and time. The 
leading-order grid-scale errors are diffusive. Eor most simulations 
reported here, we include explicit diffusion. A sufficiently large 
explicit diffusion can dominate grid-scale errors and allow the 
simulation to remain close to the true solution. However, higher- 
order grid-scale errors can introduce non-diffusive effects, such as 
dispersion. If higher-order errors project onto unstable modes, they 
can cause large differences in the solution, despite being higher 
order. The grid-scale errors in Athena respect the reflect-and-shift 
symmetry of our problem up to floating point accuracy, so even non- 
converged simulations can maintain the initial symmetry of the flow. 
In practice, we find all simulations maintain the initial symmetry, 
except simulations with Ap/po = 1 without explicit diffusion. Since 
Athena’s algorithm manifestly preserves this symmetry, we expect 
the error results from chaotic amplification of floating-point errors. 

Dedalus^ is a pseudo-spectral code (Burns et al. 2016). All 
field variables are represented as Fourier series, and the simulation 
solves for the evolution of the spectral-expansion coefficients in 
time. The code evaluates nonlinear terms on a grid with a factor 3/2 
more points than Fourier coefficients; i.e., the 2/3 de-aliasing rule. 
Lecoanet et al. (2014) (appendix D.l) describes our implementation 
of the Navier-Stokes equations. Our implementation of the dye 
evolution equation is 

^^C-Vdye [d^C + dzCz) = 

- ud^c - wCz + Vdye (dx^'^xC + d^T'c ^), (10a) 

c,-5,c = 0, (10b) 

where we use the same notation as Lecoanet et al. (2014). For 
timestepping, we use a third-order, four-stage DIRK/ERK method 
(RK443 of Ascher et al. 1997) with a total CFL safety factor of 0.6 
(i.e., 0.15 per stage). This formulation allows implicit timestepping 
of sound waves. Thus, our timestep size only adjusts with the flow 
velocity, not the sound speed. The excellent agreement between 
the highest resolution Dedalus and Athena simulations shows that 
high-wavenumber sound waves have negligible influence on the 
solution. 

The pseudo-spectral method produces almost no numerical 
diffusion. Stability concerns require explieit diffusion in nonlinear 
ealeulations. In marginally resolved simulations, diseretization er¬ 
rors manifest as Gibbs’ ringing, which is prominently visible in 
snapshots. The numerical method does not explicitly preserve the 
reflect-and-shift symmetry—numerical errors can put power into 
the asymmetric modes. However, we find that in resolved simu¬ 
lations these asymmetric modes never grow to large amplitudes. 
Thus, maintaining this symmetry gives a test for a simulation’s 
fidelity. 


^ Dedalus is available athttp://dedalus-project.org. 
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3 RESULTS 

This section describes the nonlinear evolution of the KH instability, 
provides reference solutions, and compares the performance of 
Dedalus and Athena. Section 3.1 considers unstratified simulations 
with constant initial density; both codes handle this problem eas¬ 
ily. Section 3.2 concerns simulations with a density jump across 
the shear interface. This problem shows rich behavior and poses 
significant numerical challenges. 

3.1 Unstratified simulations (Ap/po = 0) 

In this section, we discuss simulations with constant initial density 
(Ap/po = 0). Figure 1 visualizes the flow with the dye concentration 
field of the lower half of the domain for simulations with explicit 
diffusion at different resolutions and Reynolds number. Re. The 
snapshots show the state dlt - 6. Strong nonlinearity begins at 
f ~ 2, so this corresponds to at least four turnover times after the 
initial saturation of the instability. The simulations are labeled by 
the code used (A for Athena; D for Dedalus), and their horizontal 
resolution. 

The flow consists of coherent filaments of unmixed fluid with 
dye concentration close to zero or one. The filaments twist around 
the central vortex until they become thin enough to diffuse away. 
The central vortex stays coherent in all simulations, and exhibits a 
more gradual dye-concentration gradient than in the filaments. This 
reflects the smooth velocity and dye initial condition. 

3.1.1 Re = 10^ 

Many of the simulations with the same Re but different resolution 
look similar by eye. To more quantitatively assess convergence, we 
calculate the L 2 norm of the differences between dye concentration 
fields in different simulations: 



where cx and cy represent the dye concentration fields in two 
simulations, X and Y. The Athena and Dedalus grids are different, 
so we use spectrally accurate techniques to interpolate Dedalus 
solutions to the Athena grid for direct comparison (Appendix A). 
We argue in Appendix B that all simulations converge to our highest- 
resolution Dedalus simulations; thus, we assume these simulations 
are a good approximation to the “true” solution. 

Figure 2 shows the L 2 norm of the difference between dye 
concentration fields of D2048 and other simulations with Re = 10^. 
Because we believe D2048 closely represents the true solution (Ap¬ 
pendix B), we call this the L 2 norm of the error. Solutions from both 
codes approach D2048 as resolution increases. At late times, A2048 
and D1024 have roughly eight-times smaller errors than A1024 
and D512, respectively. That is, both codes exhibit third-order con¬ 
vergence. This indicates that interpolation produces the dominant 
error in Athena, which is the only third-order part of the algorithm. 
The Dedalus simulations are spatially resolved, so timestepping 
produces the dominant error source in the Dedalus simulations, 
which is also third order. We also plot errors from D512dt, which is 
run with a horizontal resolution of 512, but with half the CFL safety 
factor. D512dt is almost as accurate as D1024, showing that the 
higher accuracy of D1024 is mostly due to taking smaller timesteps. 
There are certain times (most notably near t -3.5) where the flow 
develops smaller structures, and extra spatial resolution is required. 


The errors in quantities other than dye concentration (e.g., density) 
follow similar behavior to that shown in Figure 2. 

We calculate the volume-integrated dye entropy for each sim¬ 
ulation (equation 5). Figure 3 plots the entropy as a function of 
time. Because all simulations are well resolved, there are no visible 
differences in the entropy between the different simulations. 

3.1.2 Re = 10^ 

The unmixed filaments are much thinner for Re = 10^ than for 
Re = 10^, challenging the codes. Unlike the Re = 10^ case, some 
minor visible differences appear between the solutions for Re = 10^. 
The lower-resolution simulations do not fully resolve the flow (one 
such feature is highlighted in Figure 7). 

To assess convergence, we again plot the L 2 norm of the error 
in dye concentration with respect to D2048 (Figure 4). A1024 has 
the largest errors of any simulation. At late times, the errors interact 
nonlinearly, whereas the errors in the higher-resolution Athena 
simulations stay linear and the temporal variation of the error is the 
same independent of the magnitude of the error. The ratio of errors 
of the two higher-resolution Athena simulations is about 6—in 
between second- and third-order convergence. This suggests that 
the size of interpolation errors roughly match the size of other errors 
in the code (e.g., from the Riemann problem or timestepping). 

The difference in errors between D512 and D1024 is about 
100—much larger than the difference in errors between the Athena 
simulations. D512 (not shown in Figure 1) under resolves the flow 
and includes some low-amplitude Gibbs’ ringing. Increasing the 
resolution from 512 to 1024 eliminates spatial errors because of 
the exponential convergence of spectral methods. This allows for 
very large error reduction with only modest resolution changes. The 
exponential nature of spectral methods makes convergence prac¬ 
tically binary: simulations with Gibbs’ ringing are not converged; 
simulations without Gibbs’ ringing very likely are converged. 

We plot volume-integrated dye entropy for Re = 10^ in Fig¬ 
ure 5. Like for Re = 10^, all well-resolved simulations produce sim¬ 
ilar entropy. However, the under-resolved A1024 produces slightly 
more entropy. This agrees with the heuristic that extra numerical 
diffusion leads to excess entropy generation. 

3.1.3 An effective Reynolds number ? 

We now describe Athena simulations without any explicit diffusion. 
An important question is, does the numerical diffusion in Athena 
act like an explicit diffusion? Put another way, does Athena have an 
effective Reynolds number at a given resolution for this problem? 
As we describe below and in section 3.2, the answer to this question 
is very problem dependent. 

To test this, we plot the converged volume-integrated dye 
entropy for several Reynolds numbers, along with the volume- 
integrated dye entropy for Athena simulations without explicit 
diffusion (Figure 6). The entropy evolution of N1024 is similar to 
the entropy evolution for Re = 10^. This might lead one to think 
that the effective Reynolds number of this Athena simulation is 
about 10^. 

However, a closer investigation shows that N1024 and the 
Re = 10^ simulation have different dye concentration fields which, 
by chance, result in similar volume-integrated entropies (Figure 7). 
Instead, the dye concentration field of N1024 looks like the dye 
concentration field of the (under resolved) A1024 simulation with 
Re = 10^. Figure 5 shows A1024 has a higher entropy than the 
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x/L^ x/L^ xjL^ xjL^ 


Figure 1. Snapshots of the dye concentration field in several simulations with = 0 at f = 6. The upper (lower) row shows simulations with Re = 10^ 

(10^). All the simulations with Re = 10^ are well resolved. Small differences exist between the lower-resolution Athena simulations at Re = 10^ and the 
highest-resolution Athena simulation & Dedalus simulations (e.g., near (x,z) = (0.9,0.6), see Figure 7). 


Re =10® , Ap/po=0 



Figure 2. L 2 norm of dye-concentration errors for Aplpo = 0 and Re = 10^. 
We take D2048 as the “true” solution (see Appendix B). Both Dedalus and 
Athena exhibit third-order convergence. D512dt is run with half the timestep 
size as D512. Its error is similar to D1024, showing that the higher accuracy 
of D1024 is mostly due to a smaller timestep size rather than higher spatial 
resolution. 



Figure 3. Volume-integrated dye entropy (equation 5) as a function of time 
for the four simulations with Re = 10^ shown in Figure 1. All simulations 
are well resolved, so the dye entropies are almost equal. 


true Re = 10^ solution. By removing the explicit diffusion, the 
flow evolution remains similar to A1024 (and different from the 
resolved Re = 10^ solution), but the interfaces between filaments 
are sharper, which decreases the entropy. The effects of having 
the incorrect flow field (increasing entropy), but sharper interfaces 
between filaments (decreasing entropy) happen to cancel out, so 
the entropy of N1024 is similar to that of Re = 10^. 

Although we have highlighted the differences between N1024 
and the converged solutions with Re = 10^, it is worth reiterating 
that the two solutions are in fact remarkably similar. This shows that 
N1024 roughly has an effective Reynolds number of 10^. In detail, 
however, the remaining modest differences between N1024 and the 
Re = 10^ solution demonstrate that the numerical dissipation in 
Athena is not exactly equivalent to physical dissipation via viscosity 
and thermal conduction. 


One difficulty with the notion of an effective Reynolds number 
is that it is extremely problem dependent, even at fixed resolution. In 
the next section, we introduce a small (by astrophysical standards) 
density jump into the initial condition. This completely changes 
the problem by introducing secondary instabilities which enhance 
mixing, producing very clear differences between resolved simula¬ 
tions and Athena simulations without explicit diffusion (Figure 15). 
For the constant-density problem described here, omitting diffusion 
produces less entropy. Including a density jump reverses this trend: 
simulations with only numerical diffusion undergo more mixing 
than simulations with explicit diffusion. Although assigning an 
effective Reynolds number to Athena simulations without explicit 
diffusion may be reasonably accurate for the constant-initial-density 
problem, this does not carry over to the problem with an initial den¬ 
sity jump. 
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Figure 4. L 2 norm of dye-concentration errors for Aplpo = 0 and Re = 10^. 
A1024 is not well resolved so its errors follow a different pattern than the 
other Athena simulations. The errors in A4096 are smaller than the errors in 
A2048 by ^ 6. The errors in D1024 are smaller than the errors in D512 by 
about 100. This demonstrates the fast (exponential) convergence of spectral 
methods. 
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Figure 5. Volume-integrated dye entropy (equation 5) as a function of time 
for the five simulations with Re = 10^ shown in Figure 1. The entropy of 
all simulations are very similar except for A1024; this is another indication 
that A1024 is not well resolved. 


3.2 Simulations with a density jump (Ap/po = 1) 

Both the qualitative features of the flow and the convergence prop¬ 
erties of the simulations change dramatically once we introduce 
an initial density jump (Ap/po ^ 0). Unlike the unstratified case, 
secondary instabilities of the filaments produce small-scale struc¬ 
tures in the flow. These secondary instabilities, and the resulting 
small-scale features, depend on the resolution and the code used. 
As a result, simulations with a nonzero density jump require far 
more computational resources than the unstratified simulations pre¬ 
sented in the previous section. We limit the simulations with explicit 
diffusion to Re = 10^—our finite-computing budget precludes solu¬ 
tions for Re = 10^. The largest simulations required roughly 10^ 
core-hours. 

Figure 8 shows the dye concentration for different simulations 
at different times. In both Dedalus simulations, and the highest- 
resolution Athena simulation, the outer filaments (i.e., those outside 
the central vortex) become unstable to a sausage-like mode (see the 
panel in Figure 10 for an example). Lower-resolution Athena simu¬ 
lations also undergo a separate instability of the inner filaments of 
the vortex. We refer to these two instabilities at the outer-filament in¬ 
stability (OFI) and the inner-vortex instability (IVI) (see Figure 10 


Figure 6. Volume-integrated dye entropy (see section 2.2) as a function of 
time with Ap/po = 0, for three resolved simulations with different Re, as 
well as three Athena simulations with no explicit diffusion (dashed lines; 
labeled with N, for no explicit diffusion, and their horizontal resolution). The 
entropy of N1024 and the simulation with Re = 10^ are very si mil ar. Their 
flow fields show minor differences (see Figure 7). Note that the entropy 
decreases with increasing resolution in the simulations without explicit 
diffusion. This is not the case in simulations with an initial density jump 
(see Figure 16). 



0.90 0.92 0.94 0.90 0.92 0.94 

xjL^ xlL^ 


Figure 7. Snapshots of the dye concentration field between 0.89 < x < 
0.95 and 0.55 < z < 0.61, at r = 6 for Ap/po = 0. All simulations 
use Athena, either with Re = 10^ (left column) or no explicit diffusion 
(right column). The three rows have different resolutions. This zoom-in 
of Figure 1 highlights the differences between simulations at different 
resolutions—however, for the most part, the simulations look very similar. 
A2048 & A4096 represent resolved simulations with Re = 10^. Although 
the entropies for N1024 (upper right plot) & A4096 (lower left plot) track 
each other (Figure 6), the dye concentration fields exhibit minor differences. 
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x/L^ xjL^ xjL^ x/L^ xjL^ 


Figure 8. Snapshots of the dye concentration field in several simulations with Aplpo = 1 and Re = 10^. Each row corresponds to a different time. The 
low-resolution Athena simulations suffer from a secondary instability (seen at f = 4) in the middle of the vortex, which is not present in the Dedalus simulations 
nor A16384. This causes substantial differences at later times. A16384 and both Dedalus simulations stay very similar at late times, although small differences 
develop from chaos (see section 3.2.2). 


for examples). These instabilities are similar to the baroclinic sec¬ 
ondary instabilities discussed in Reinaud et al. (2000); Fontane & 
Joly (2008). The competition between these two instabilities plays 
a crucial role in the evolution of the system. 

We plot the L 2 norm of the error in dye concentration with 
respect to D4096 in Figure 9. As described in Appendix B, we be¬ 
lieve D4096 approximates the true solution. The difference between 
D3072 and D4096 are smaller than the differences between any 
other pair of simulations. At later times, even the errors between 
D3072 and D4096 become large. In section 3.2.2 we attribute this 
late-time behavior to chaos. 

Figure 9 shows that at early times, the low-resolution Athena 
simulations diverge exponentially from D4096 with an inferred 
growth rate of about 8. The IVI produces this divergence. Fur¬ 
thermore, the four Athena simulations with resolutions between 
1024 and 8192 are all equally spaced horizontally in Figure 9. The 
horizontal-axis spacing is log 2/2 time units. This suggests that the 
same instability exists independent of resolution, but the amplitude 
of the perturbation that seeds the instability drops by 16 when the 
resolution doubles. Though numerical errors seed the growth, the 


constant growth rate of the IVI suggests it is a physical instability 
(we demonstrate this in section 3.2.1). 

The IVI is a robust feature of low-resolution Athena simu¬ 
lations. Using the Roe integrator, second-order reconstruction, or 
shifting the initial condition by half a grid point does not affect 
the development of this instability (as confirmed using the L 2 er¬ 
ror), but can cause visible differences in the flow evolution. This 
demonstrates that grid-scale errors drive the IVI. Using first-order 
reconstruction suppresses the IVI, but the enhanced numerical dif¬ 
fusion causes large errors. We have also tried adding low-amplitude 
(up to 10""^) white noise to the initial density or pressure. These do 
not cause any visible changes to the IVI. The flow forgets some of 
the detailed information of its initial condition (see section 3.2.3). 

The highest-resolution Athena simulation (A16384) does not 
develop the IVI. This demonstrates that the initial condition is in 
fact stable to the IVI; the problem is well-posed. Rather, numerical 
errors seed the IVI at some later time, during the evolution of the 
flow. Although some numerical errors are still inevitably present, 
A16384 does not develop the IVI because the “base state” of spi¬ 
ralling filaments of unmixed fluid also succombs to the OFF In this 
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Figure 9. L 2 norm of dye-concentration errors for Aplpo = 1 and Re = 10^. 
D3072 and D4096 are the closest pair of simulations, suggesting that D4096 
is a good approximation to the true solution. All Athena simulations except 
A163 84 diverge away from D4096 exponentially with a rate of 8, suggesting 
the growth rate of the inner vortex instability (see Figure 10) is also 8. 
The errors in the lower-resolution Dedalus simulation and A163 84 grow 
exponentially with a rate of about 2-3. We interpret this divergence as due 
to chaos (see section 3.2.2). D3072 has errors smaller than D2048 by ^ 4, 
consistent with third-order convergence set by our choice of timestepping 
algorithm. 


case, the OFI disrupts the inner vortex before the IVI grows to large 
amplitudes (see Figure 10). 

The absence of the IVI is a robust feature of our Dedalus sim¬ 
ulations. We confirmed the stability of the base state by re-running 
D2048 with low-amplitude white noise added to the initial condi¬ 
tion; we also re-initialized D2048 from the Athena initial condition. 
This introduces small but non-random grid-representation differ¬ 
ences (section 3.2.1). In both cases, we recover the same evolution. 
However, we can trigger the IVI in Dedalus with a large (~ 10% by 
energy) perturbation to the initial condition (section 3.2.3). 

Figure 10 summarizes the relation between the two seeondary 
instabilities in this problem. For a constant initial density (left 
panel), the system evolves toward a stable state characterized by 
spiraling filaments. Small differences in initial conditions, integra¬ 
tion algorithms, presence of dissipation, etc., cause only minor 
changes in the evolution. We hypothesize that a similar spiral state 
also exists for Ap/po = 1, and that it could be reached from some 
initial condition IC’. However, our simulations demonstrate that 
the spiral state is now unstable. Thus, small errors lead to the large 
differences in evolution. 

Small perturbations to the hypothetical IC’ of Figure 10 would 
lead to trajectories that either develop the OFI or the IVI. However, 
our chosen initial condition, IC, is squarely in the attracting basin 
of the OFI. Thus, infinitesimal perturbations to IC will still lead 
to the OFI. Errors introduced by numerical hydrodynamics cause 
the codes to not follow the correct trajectory (solid black line). 
Certain types of errors can cause trajectories to diverge from the 
correct solution, sometimes toward the IVI (dashed grey lines). 
Alternatively, sufficiently large initial perturbations can also knock 
the system into the attracting basin of the IVI (section 3.2.3). 

We note that the phase space for this problem is very high 
dimension, and that the outer filament instability and inner vortex 
instability represent two (likely non-parallel) unstable directions of 
the spiral state’s stable manifold. Thus, both instabilities can act 
simultaneously, which sometimes occurs in simulations. 

Figure 11 shows the volume-integrated dye entropy of the sim¬ 
ulations shown in Figure 8. The entropy follows a similar evolution 


in every simulation. To visualize the small deviations, the bottom 
panel shows the entropy with reference solution D4096 subtracted 
off. All the simulations diverge from D4096, but more accurate 
simulation diverge later, with D2048 and A16384 developing small 
differences later than any other simulation. The relation between 
entropy and resolution is more complicated for Ap/po = 1 than for 
Ap/po = 0 (Figures 3 & 5). 

Apart from the dye concentration field, many of the other flow 
quantities follow similar patterns. Figure 12 shows several quanti¬ 
ties from D4096 at ^ = 6. The mass density is almost the inverse 
of the dye concentration. This indicates that compression is not an 
important part of the large-scale dynamics. Lacking mass diffusion, 
the density shows sharper gradients than the concentration field. 
Temperature diffusion and rapid sound waves regularize the density 
evolution. These effeets limit large temperature gradients, and keep 
the flow in local pressure equilibrium. 

The velocity divergence field is characterized by a large scale 
quadrupole centered at the vortex, and large amplitude, small scale 
features near the boundaries of filaments. The most prominent 
feature of the vorticity field is the central vortex, which is a remnant 
of the initial shear. Small-seale vortex sheets and filaments perhaps 
result from the incomplete roll-up of the initial condition due to 
secondary instabilities. 

Throughout this paper, we compare different solutions by cal¬ 
culating the L 2 norm of the difference between dye concentration 
fields. We have made similar comparisons between simulations 
with Re = 10^ and Ap/po = 1 using the Li norm of the difference 
between dye concentration fields, and using the L 2 norm of the dif¬ 
ference between the three other fields shown in Figure 12. We find 
the results to be qualitatively similar in all cases. This is expected 
given the similarity between the fields. 


3.2.1 Inner-vortex instability 

To determine the origin (physical vs numerical) of IVI, we initialize 
a Dedalus simulation with horizontal resolution 2048 with the out¬ 
put from A2048 at ? = 3.2. We call this simulation D2048r. Figure 9 
shows that A2048 is still in the linear phase of the IVI at this time. 
In Figure 13, we plot the dye eoncentration field at ? = 3.2 and t = 4 
for D2048, A2048, and D2048r. At ? = 3.2, the simulations all look 
the same. However, the instability becomes nonlinear by ^ = 4, pro¬ 
ducing large changes in the dye concentration field. D2048 shows 
no signs of the IVI. However, D2048r looks almost identical to 
A2048. The L 2 norm of the difference of dye concentration fields 
between D2048r and D4096 almost exactly follows the norm of the 
difference between A2048 and D4096. 

This shows that the IVI is a physieal instability of this system. 
It is not seen in the Dedalus simulations or the highest-resolution 
Athena simulation because the initial eondition does not projeet suf¬ 
ficiently onto its unstable modes. Errors in low resolution Athena 
simulations incorrectly excite perturbations unstable to the IVI. 
Dedalus simulations, and the highest-resolution Athena simulation, 
suppress noise well enough the instability never becomes nonlin¬ 
ear. In our phase-space diagram (Figure 10), the lower-resolution 
Athena simulations do not properly follow the black line, and in¬ 
stead meander to the right, becoming unstable to the IVI. D2048r is 
initialized to the right of IC’, so it develops the IVI just like A2048. 

As a final test, we started a Dedalus simulation from the output 
of an Athena simulation at ^ = 0. This tests whether dynamical 
evolution causes the IVI, rather than small differences between the 
implementation of the initial conditions. Although this introduced 
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Figure 10. Schematic phase-space diagram for Ap/po = 0 (left) and Ap/po = 1 (right). For constant initial density, the system has a stable state with 
ever-narrowing spiral filaments. We hypothesize that there is an initial condition IC’ (right panel) leading to a similar spiral state for Ap/po = 1. But this state 
is now unstable to the outer filament instability (OFI) and the inner vortex instability (IVI). Our chosen initial condition’s (IC) trajectory (solid black line) 
approaches the spiral state, but becomes unstable to the OFI. Errors introduced by the numerical hydrodynamics may cause deviations in the trajectory leading 
to the IVI (dashed grey lines). 
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Figure 11. Volume-integrated dye entropy (equation 5) as a function of 
time for simulations with Aplpo = 1 and Re = 10^. The top panel plots the 
entropy, and the bottom panel plots the entropy deviation from D4096. The 
entropy of all the simulations diverge from D4096, but the less-accurate 
simulations diverge faster. For each Athena simulation, the entropy initially 
increases faster than D4096 when it starts to diverge. At later times, the 
entropy sometimes drops below the entropy of D4096. 


root mean squared differences in the horizontal velocity of ^ 4 X 
10""^ at ? = 0, the Dedalus simulation did not develop the IVI. 


3.2.2 Chaos 

At around ^ ^ 4, D2048, D3072, and A16384 start to diverge 
exponentially from D4096 (Figure 9). The differences increase with 
a growth rate of about 2-3, much lower than the growth rate of 8 
of the IVI found in the lower-resolution Athena simulations. We 
interpret the differences between the simulations as due to chaos. 
The faster divergence discussed in section 3.2.1 is inconsistent 
with chaos since it is resolution dependent and only seen in low- 
resolution Athena simulations. 

A system is chaotic if small differences between initial condi¬ 
tions grow exponentially in time. To confirm the system is chaotic, 
we calculate a ‘Tocal-in-time” Lyapunov exponent (i.e., growth 
rate). We pick a time and simulation, and look for linearly unstable 
perturbations. This requires solving an eigenvalue problem. The 
largest unstable eigenvalue is the Lyapunov exponent. Appendix C 
details this procedure. 

This calculation does not inculde base-state time evolution 
(i.e. we consider a “local-in-time” calculation). The most unstable 
eigenvector at a time might differ significantly from the most 
unstable eigenvector at a nearby time + At. Then it would be 
impossible for perturbations to grow at the Lyapunov exponent over 
times ~ At. We interpret our “local-in-time” Lyapunov exponents 
as an upper bound on the growth rate of perturbations due to chaos 
(up to logarithmic corrections), and as a heuristic measure of the 
strength of chaos in this problem. 

We calculated the Lyapunov exponent for D2048 with Re = 
10^ and Ap/po = 1 at two times, t = 2.5 and t = 4.5. We find 
Lyapunov exponents of 3t=2.5 ~ 2.1, and 4^=4 5 ^3.7. Thus, the 
exponential growth of differences between either D2048, D3072, or 
A16384 and D4096 is consistent with chaos. However, the growth 
rate of the differences between the lower-resolution Athena simula¬ 
tions and D4096 is much larger than the Lyapunov exponent. These 
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Figure 12. Plots of dye concentration (c), mass density (p), the divergence 
of the velocity (V • u), and the vorticity (co = • V x m) in D4096 with 

Aplpo = 1, Re = 10^ at r = 6. The divergence of the velocity and the 
vorticity are measured in units of Wflow/l-x- The dye concentration and mass 
density fields are almost inverses of each other. The divergence of the 
velocity is largest at the interfaces between filaments, whereas the vorticity 
shows the location of vortices. 



Figure 13. Snapshots of dye concentration field for Re = 10^ and Ap/po = 
1. D2048r is a Dedalus simulation restarted with the A2048 output att = 3.2. 
At this time, the inner vortex instability is still in the linear phase, so there 
are no visible differences between the three simulations. Att = 4, the IVI 
is very nonlinear, producing large differences between D2048 and A2048. 
This instability also takes place in D2048r, and the dye concentration fields 
of A2048 and D2048r are nearly identical. This demonstrates that the IVI is 
physical, but is seeded by errors in the lower-resolution Athena simulations 
that are not present in the Dedalus simulations or the highest-resolution 
Athena simulations. 



Figure 14. Snapshots of dye concentration field for Re = 10^ and 
Ap/po = 1- D2048p is a Dedalus simulation with an initial vertical ve¬ 
locity that includes power over a range of Fourier modes (equation 12), in 
contrast to the single mode initial conditions focused on throughout the 
rest of this paper. At f = 2 all solutions look the same, indicating that 
longest wavelength mode has the largest growth rate. At f = 4, D2048p has 
developed the IVI, as well as other deviations from the Dedalus & Athena 
simulations away from the vortex. 


differences are inconsistent with chaos, instead being due to the IVI 
(section 3.2.1). 

The simulations with Ap/po = 0 do not appear to diverge 
from one another in the same way. The highest-resolution Dedalus 
simulations converge at late times. We also calculate the Lyapunov 
exponent for D1024 with Re = 10^ and Ap/po = 0 at t = 6. We 
find At=^ ^ 0.4. Although this seems inconsistent with our finding 
that the Dedalus simulations approach each other with time, recall 
that this ‘Tocal-in-time” calculation gives an upper bound on the 
growth rate due to chaos (up to logarithmic corrections). Because 
the turnover time is 1, a Lyapunov exponent less than 1 suggests 
that small perturbations cannot grow before the background state 
changes substantially. To show definitively that the Ap/po = 0 
solution is not chaotic, one should maximize the amplification of an 
initial perturbation over several turnover times, e.g., between t = 6 
and t = 9 (for instance, using the adjoint method, e.g., Kerswell 
et al. 2014). 


3.2.3 Initial condition 


Although our chosen initial condition does not lead to the IVI for 
converged simulations, one might wonder if other initial conditions 
do lead to this instability. We performed several Dedalus simula¬ 
tions that add low-amplitude white noise to the initial condition 
(e.g., see section 3.2.1). None of these simulations develop the IVI. 

We now consider a simulation in which we include perturba¬ 
tions to the initial condition with order unity amplitude and large 
wavelengths. Equations 8 still hold for all quantities except the 
vertical velocity, which we now take to be 


Uz- A (sin(2:7rx) -f f{x)) X 


exp 


(z-Zl)^ \ 
0-2 ] 


+ exp 


(z - Zlf \ 

0-2 jj’ 
( 12 ) 


where f{x) includes Fourier modes two-ten. Each mode receives a 
random phase and random amplitude uniformly distributed between 
-0.05 and 0.05. Thus, f(x) represents about a 10% perturbation to 
the single sine mode initial condition. 

Figure 14 shows snapshots of the dye concentration field for 
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Figure 15. Snapshots of dye concentration field for Ap/po = 1. N4096 is 
an Athena simulation with no explicit diffusion. For comparison, we also 
plot D4096 (Re = 10^). Secondary instabilities occur very early at many 
locations in N4096. By f = 6, the simulation has broken its initial symmetry 
(we only plot the bottom half). The secondary instabilities produce signifi¬ 
cant mixing, leading to greater entropy generation than in simulations with 
explicit diffusion (Figure 16). 

this simulation, denoted D2048p, along with D2048 and A2048 
for comparison. Kit -2, all three simulations look identical. This 
indicates that the lowest wavenumber Fourier mode grows faster 
than the other modes included in our initial condition. 

By r = 4 the perturbations from the other Fourier modes 
produce significant changes to the dye concentration field. D2048p 
now displays the IVL In addition, large differences appear away 
from the vortex, where the Dedalus and Athena simulations look 
almost identical. Because the new initial condition does not respect 
the shift-and-reflect symmetry of the problem, the two half domains 
have different features (we only show the bottom half). 

3.2.4 Simulations without explicit diffusion 

Lastly, Figure 15 compares the resolved simulations at Re = 10^ 
with an Athena simulation with horizontal resolution 4096 without 
explicit diffusion (N4096). The simulation without explicit diffu¬ 
sion exhibits many secondary instabilities early in the evolution 
(between t - 2 and t = 4). Unlike the lower-resolution simulations 



Figure 16. Volume-integrated dye entropy (equation 5) as a function of 
time for simulations with Ap/po = 1. D4096 is run at Re = 10^, and all 
simulations labeled with N are run with Athena with no explicit diffusion. 
At early times, the highest-resolution runs without explicit diffusion have 
the lowest entropy. However, at around f = 5, the lower-resolution runs 
without explicit diffusion have lower entropy. D4096 has the lowest entropy 
at late times. This indicates that simulations without explicit diffusion have 
greater numerical mixing compared to simulations with explicit diffusion. 
This becomes more prominent as the resolution increases. By contrast, in 
the simulations without an initial density jump, explicit diffusion leads to 
more mixing, and for simulations without explicit diffusion, increasing 
resolution decreases mixing (Figure 6). 

at Re = 10^, the secondary instability is not limited to the IVL 
Instead, instabilities grow throughout the domain at locations of 
strong shear. 

These instabilities shred apart the vortex, leading to vigorous 
mixing. Figure 16 compares the volume-integrated dye entropy of 
Athena simulations with no explicit diffusion at different resolutions 
with D4096. Simulations without explicit diffusion produce almost 
no entropy until f ^ 3.5. At this time, the secondary instabilities 
start to cause diffusion at the grid-scale. This generates entropy 
more rapidly than the explicit diffusion of D4096 (or any of the 
other simulations with explicit diffusion). For ? > 5, the entropy of 
the simulations without explicit diffusion is larger than the entropy 
of D4096. Paradoxically, the entropy increases as the resolution 
increases. Our expectation is that the entropy generation should 
decrease as Re increases. However, we do not have any resolved 
simulations with higher Re for comparison, so we cannot present 
evidence that this additional mixing is spurious. But this problem 
shows that introducing an explicit diffusion in Athena can decrease 
the diffusion in the simulation. 


4 CONCLUSION 

This paper describes several converged, nonlinear solutions to the 
Kelvin-Helmholtz (KH) problem. By using a smooth initial condi¬ 
tion and explicit diffusion, we demonstrate that solutions remain 
virtually identical (for constant initial density) or very similar (for 
an initial density jump of one) with resolution above a certain 
threshold. This permits a well-defined reference solution for this 
problem, against which errors can be accurately estimated. We ver¬ 
ify this using two codes, Dedalus and Athena, with very different 
numerical methods (pseudo-spectral and Godunov, respectively). 
Previous KH test problems either did not use smooth initial condi¬ 
tions, or did not include explicit diffusion. Absent these two choices, 
the KH problem cannot be quantitatively compared between codes 


© 2015 RAS, MNRAS 000, 1-15 














12 Lecoanet et al 


because the solutions depend sensitively on grid-scale errors and 
do not converge with increasing resolution. 

We first study simulations with a constant initial density (sec¬ 
tion 3.1). We find converged solutions to this relatively easy prob¬ 
lem with Reynolds numbers (Re) as high as 10^. The solution is 
characterized by the continual roll-up of the initial vortex sheet, 
producing alternating filaments of unmixed material (Figure 1). 
We find third-order convergence in both Dedalus & Athena for 
simulations with Re = 10^ (Figure 2), and better than second- 
order convergence in both codes for simulations with Re = 10^ 
(Figure 4). 

To quantify mixing in the simulations, we calculate the 
volume-integrated dye entropy as a function of time for several 
Reynolds numbers, as well as for Athena simulations without ex¬ 
plicit diffusion (Figure 6). As the Reynolds number increases, the 
entropy generation decreases monotonically. Similarly, as the reso¬ 
lution of Athena simulations without explicit diffusion increases, 
the entropy generation also decreases monotonically. The entropy 
of one Athena simulation without explicit diffusion is very close 
to the entropy of the Re = 10^ simulation, although the solutions 
show minor differences (Figure 7). These small differences indicate 
that the numerical diffusion in Athena does not act precisely as 
a physical diffusion from viscosity and/or thermal conductivity. 
For certain applications however, assigning an effective Reynolds 
number to ideal fluid simulations may suffice. This does not appear 
to be the case for KH simulations with density jumps, as we now 
discuss. 

Including an initial density gradient aligned with the velocity 
gradient makes the problem much richer (section 3.2). The rolled-up 
vortex-sheet filaments becomes unstable in at least two ways: the in¬ 
ner vortex instability, and/or the outer filament instability (Figures 8 
& 10). The Dedalus simulations and highest-resolution Athena 
simulation only exhibit the outer filament instability, whereas the 
lower-resolution Athena simulations also exhibit the inner vortex 
instability. Adding small amplitude noise to the initial condition 
does not produce the inner vortex instability in Dedalus, demon¬ 
strating that our chosen initial condition is not susceptible to this 
instability; instead, numerical errors seed the inner vortex instability 
throughout the evolution of the Athena simulations. It is not sur¬ 
prising that Dedalus is more accurate than Athena for this smooth 
flow—the Godunov method is designed for simulating flows with 
shocks. However, it is not well appreciated that the pseudo-spectral 
method is able to solve the full Navier-Stokes equations with Mach 
number order unity. 

We use the L 2 norm to quantify the difference between dye 
concentration fields of different simulations, and find the inner 
vortex instability grows at a rate of ^ 8, independent of resolu¬ 
tion (Figure 9). Furthermore, a Dedalus simulation initialized with 
an Athena state in the linear phase of the inner vortex instability 
develops the instability in the same way as Athena (Figure 13), 
demonstrating the physical, rather than numerical, nature of the 
instability. 

Adding a large (~ 10% by energy) perturbation with multiple 
Fourier modes to the initial velocity in Dedalus can seed the inner 
vortex instability (section 3.2.3). Although this suggests that the 
inner vortex instability is possibly generic for KH instabilities in 
astrophysics, we believe the single-mode initial condition discussed 
throughout the rest of this paper is still particularly valuable for 
a test problem. Because small numerical errors can produce large 
differences in the solution, one can assess by eye the fidelity with 
which a code is solving the fluid equations. This KH test problem 


is difficult, which we believe makes it interesting. In contrast, an 
unresolved KH problem is not a good test of fluid codes, because 
noise due to numerical errors can masquerade as higher-fidelity 
solutions. 

The Dedalus simulations and highest-resolution Athena sim¬ 
ulation also diverge from each other exponentially at late times, 
but with a much smaller growth rate ^ 2 - 3. In section 3.2.2 we 
calculate the maximum Lyapunov exponent of the flow, and argue 
that chaos drives the divergence. The Lyapunov exponent represents 
the maximum possible rate of divergence of solutions due to chaos 
(up to logarithmic corrections). At late times when the Dedalus 
simulations and highest-resolution Athena simulation begin to di¬ 
verge, the Lyapunov exponent is ^ 3.7, so the divergence we see is 
consistent with chaos. Because the system is chaotic, our solutions 
are not as accurate as the solutions with constant initial density. 
We still find power-law convergence in the Dedalus simulations at 
fixed time (Figure 9). However, the amount of time that a solution 
maintains a fixed level of accuracy increases only logarithmically 
with resolution. 

For the initial condition with a density jump, we also compare 
a high-resolution Athena simulation without explicit diffusion to 
our converged (within the limits of chaos) simulations with Re = 
10^. Secondary instabilities pervade the simulation without explicit 
diffusion (Figure 15). The secondary instabilities cause enhanced 
mixing, and at late times, the simulations without explicit diffusion 
have higher entropy than the Re = 10^ simulation (Figure 11). 
Introducing explicit diffusion into Athena can reduce the diffusion 
in the simulation. For this reason, we hypothesize (but cannot prove) 
that this small-scale structure is likely unphysical, and would not 
develop for any reasonable initial condition or Reynolds number. 
This highlights that a solution with more small-scale structure is 
not necessarily better. 

Although we only describe simulations with an initial density 
ratio of one, we have experimented with larger initial density ratios 
(e.g., 4). Preliminary investigation suggests that vigorous secondary 
instabilities become increasingly prominent as the density ratio 
increases, greatly enhancing mixing. Though it’s a common practice 
to leave out explicit dissipation to model the high Reynolds numbers 
relevant in astrophysics, our results suggest that including explicit 
diffusion may provide a very effective way to reduce diffusion in 
astrophysical simulations with very large density ratios. We stress 
that these large density ratios are common in astrophysical problems 
such as star formation or galaxy formation. Our results demonstrate 
just how subtle and computationally challenging it is to correctly 
capture mixing in these environments (even restricting ourselves to 
hydrodynamics, which is likely a poor approximation). 

There are many remaining questions left unanswered in this 
paper. It is unclear how the Athena algorithm seeds the inner vortex 
instability. We did not search for the critical perturbation amplitude 
that will cause a Dedalus simulation to exhibit the inner vortex insta¬ 
bility. Because of limited computer time, we did not find converged 
Dedalus or Athena simulations with Ap/po = 1 and Re = 10^. 
Perhaps, contrary to expectation, increasing the Reynolds number 
of the system does increase the entropy production, as found in the 
Athena simulations without explicit diffusion. Future work should 
also test the Galilean invariance of these simulations, test initial 
conditions with an interface at an angle to the grid, and extend this 
analysis to larger density ratios. 

We hope this study provides a well-posed test problem for 
future codes used in astrophysics. It would be valuable to carry 
out this test problem with unstructured/meshless methods (e.g.. 
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Springel 2010; Duffell & MacFadyen 2011; Hopkins 2015) to un¬ 
derstand their convergence properties on this challenging problem. 
Toward this goal, we include the reference solutions to these KH 
problems in the supplementary material accompanying this paper. 
Introducing smooth initial conditions and explicit diffusion allows 
us to calculate a converged reference solution and compare between 
codes. The competing secondary instabilities for initial conditions 
with a density jump of one provides a stringent test of the fidelity 
with which a code solves the Navier-Stokes equations, making it a 
great test problem. 
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APPENDIX A: INTERPOLATION TO A COMMON GRID 

The grid points used in Dedalus and Athena differ slightly. Por a 
periodic simulation between 0 and L with spacing Av, the Dedalus 
grid points are {0, Av, 2Av,..., L - Av}, whereas the Athena grid 
points are {Av/2,3Av/5,..., L - Av/2}. We use two spectrally ac¬ 
curate methods for interpolating Dedalus and Athena data to a 
common grid. In several cases we test both methods and find excel¬ 
lent agreement. 

Our first method is spectral interpolation. The Dedalus data 
can be viewed as either N = Ljl^x values on grid points or N 
Fourier coefficients. We can pad the Fourier coefficients with zeros 
and transform to a grid of any uniform spacing. Going from N 
to 2A points, we can compare every other entry to the Athena 
data. In a second method, we multiply the Fourier coefficients by 
exp(/k;(;Av/2). A Fourier transform then shifts the grid points by 
Av/2, to align the Dedalus grid with the Athena grid. We follow 
the same procedure in the z direction. 

Throughout this paper, we treat the Athena data (including 
initial conditions) as cell-centered data. However, the data are actu¬ 
ally volume-averaged. The lowest-order differences between cell- 
centered and volume-averaged quantities scales as ~ Av^. Thus, any 
errors associated with these differences should decrease with order 
2. In all cases studied here (i.e.. Figures 2, 4, & 9), we find better- 
than second-order convergence. This suggests that differences due 


to interpreting data as cell-centered rather than volume-averaged is 
not the dominant source of error. 


APPENDIX B: CONVERGENCE TO A “TRUE” 
SOLUTION 

This paper describes a series of calculations of the nonlinear evo¬ 
lution of the KH instability, as a function of resolution and Re. 
Without an analytic solution, we must assess the quality of the 
solutions carefully. We make two assumptions to help interpret our 
results. 

(i) Dedalus and Athena converge to the same solution at fixed Re 
as the resolution increases. We refer to this unattainable “Platonic 
ideal” solution as the true solution. 

(ii) The distance (given a choice of norm) between two solu¬ 
tions at different resolutions (for the same code), is larger than the 
distance between the higher-resolution simulation and true solution. 

Our simulations support these assumptions, but it is very diffi¬ 
cult, if not impossible, to prove these statements. The existence and 
uniqueness of solutions to the Navier-Stokes equations remains an 
active field of research (Fefferman 2000). 

To support these assumptions. Figure B1 plots the relative 
differences between simulations with Re = 10^ and Ap/po = 0 at 
t = 6 (described further in section 3.1.2). The top panel assumes 
our highest-resolution Dedalus solution is the true solution. The 
bottom panel assumes our highest-resolution Athena solution is the 
true solution. To assess the deviations, we plot the norm of the 
difference of dye concentration fields. This allows us to define an 
error (alternatively a distance) between two solutions X and Y as 

KX,Y)=L2(Cx-Cy), (Bl) 

where Cx and Cy are the dye concentration fields of solutions X 
and Y, respectively. Figure B1 remains mostly unchanged if we 
compare lower Reynolds number simulations with Re = 10^ and 
^P/Po = 0^ although the picture is more complicated for Ap/po = 1 
due to chaos (see section 3.2). 

Both the Athena simulations and the lower-resolution Dedalus 
simulations are converging to D2048. The top panel of Figure Bl 
therefore suggests a true solution lives very close to D2048 (as¬ 
sumption (i)). The Athena simulations converge slower than the 
Dedalus simulations because in Dedalus spatial errors decrease 
exponentially. 

The bottom panel of Figure B1 shows that A4096 is a worse 
approximation to the true solution. This is because the Dedalus 
simulations are not converging to A4096. 

One could argue that perhaps the Athena simulations are con¬ 
verging to a solution near A4096 and the Dedalus simulations are 
converging to a different solution near D2048. However, this would 
require the error of the Athena simulations with respect to D2048 to 
stay constant as the resolution increases, contrary to the top panel. 
Thus, we believe that both codes are converging to a true solution 
close to D2048 (assumption (i)). 

Presumably if Athena were run at very high resolutions, it 
would become closer to the true solution than D2048. In this case, 
we hypothesize that both Dedalus and Athena simulations would 
converge to this very high-resolution Athena simulation. For the 
range of resolutions examined in this paper, our highest-resolution 
Dedalus simulation is always closest to the true solution. 

The main idea behind assumption (ii) is the convergence prop¬ 
erties of the algorithms used in Athena and Dedalus. Specifically, 
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Figure Bl. Differences between different solutions for Re = 10^ and 
Ap/po = 0 at f = 6. In the top panel, the Dedalus simulation with hor¬ 
izontal resolution 2048 (D2048) is assumed to be the true solution, and in 
the bottom panel the Athena simulation with horizontal resolution 4096 
(A4096) is assumed to be the true solution. The error with respect to the 
assumed true solution is the L 2 norm of the difference of the dye concen¬ 
tration fields (equation Bl). The top panel shows that both Athena and 
Dedalus are converging to the high-resolution Dedalus solution, supporting 
the assumption that it is close to the true solution. In the bottom panel, the 
Athena solutions are converging to the high-resolution Athena simulation, 
but the Dedalus solutions are not. This suggests that the Athena solution is 
further from the true solution than the Dedalus solutions. 


both codes are better than first-order accurate. Imagine we somehow 
know the true solution to our problem, T. If we run a high-resolution 
simulation, SI, we calculate the error, 

e(Sl,T) = Ei. (B2) 

Now suppose we run another simulation S2 at double resolution. If 
SI and S2 are converging to T, then 

e(S2,T) = £ 2 < y, (B3) 

where better-than first-order accuracy implies the inequality. Athena 
is between second- and third-order accurate, so we expect Ei/4 < 
E 2 < Eil^ for Athena. Dedalus is exponentially accurate in space, 
and third-order accurate in time. Thus, for Dedalus, we should 
expect E 2 < £^ 1 / 8 . Nevertheless, equation B3 implies, via the 
triangle inequality, 

e(S 1 , S2) > y > e(S2, T), (B4) 

which shows that assumption (ii) holds. One can check visually that 
equations B3 & B4 hold for the simulations described in Figure Bl, 
assuming that T is very close to D2048. 


APPENDIX C: LYAPUNOV EXPONENT CALCULATION 


One can write the equations of motion (equations 1) as 

= E{U\ (Cl) 


where U = (p, m, E) is the state vector. Then infinitesimal perturba¬ 
tions to U evolve according to the equation 




SU, 


(C2) 


where 6E16U is the Frechet derivative, evaluated at U(t). 

To calculate the “local-in-time” Lyapunov exponent, we fix 
the state vector to its value at a specific time t = to. The maximum 
Lyapunov exponent is the greatest eigenvalue of {6Ej5U)u{tQ)- It is 
impractical to solve this eigenvalue problem directly—a 2D prob¬ 
lem with resolution greater than 1000 in each direction generates 
very large matrices. Instead, we solve an initial value problem by 
picking bC/(T = 0 ), and evolving 


dr5U = 


6E 


SU, 


U(tQ) 


(C3) 


where r should not be thought of as time, as we have fixed the 
background state U iot - to. The maximal Lyapunov exponent is 


A - lim log 

T—>00 


IIAu(t)|| \ 

l|AC/(0)||/’ 


(C4) 


for some norm || • ||. We choose ^/\\u\p. This is equivalent to the 
power method. 

We solve equation C3 in Dedalus using two methods. Both 
methods give very similar Lyapunov exponents. In the first, we 
directly evolve the linearized equations C3. We treat terms inde¬ 
pendent of C/o implicitly, and treat all other terms explicitly. The 
second method uses an iteration. On each iteration, we evolve the 
full equations of motion (equation Cl) for Uo + 6Ui for a time 
At to to get a state we call Ui(to + At). The initial perturbation for 
the next iteration becomes 6Ui+i oc Ui(to + At)-AUo, but with norm 
10"^. A C/o = U(to + At) - U(to) is the change in the unperturbed 
solution Uo over the time At. We normalize after each iteration to 
ensure the perturbations stay linear. 

In both cases, we initialize the calculation with a guess of 
random noise. After substantial evolution, the system begins to 
execute limit cycles (or seems to be close to a limit cycle for 
Ap/po = 1 at C = 2.5). We report the growth averaged over a 
limit cycle. 
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